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Abstract. We investigate the dynamics in a galactic potential with two reflection 
symmetries. The phase-space structure of the real system is approximated with 
a resonant detuned normal form constructed with the method based on the Lie 
transform. Attention is focused on the stability properties of the axial periodic 
orbits that play an important role in galactic models. Using energy and ellipticity as 
parameters, we find analytical expressions of bifurcations and compare them with 
numerical results available in the literature. 
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1. Introduction 

To determine salient features of the orbital structure of non-integrable 
potentials is an important topic in dynamical astronomy. Techniques 
based on the various versions of perturbation theory have been applied 
to several examples and with various degrees of approximation (for 
a review, see, e.g., Contopoulos, 2002). Of particular interest is to 
understand motion in potentials which seem to be suitable to describe 
elliptical galaxies. Among other features, the knowledge of the stability 
properties of the main periodic orbits is of paramount importance, since 
the bulk of density distribution is shaped by the stars in regular phase- 
space regions around stable periodic orbits (Binney and Tremaine, 
1987). In particular, for triaxial ellipsoids, periodic orbits along sym- 
metry axes play a special role. An enormous effort has therefore been 
devoted to investigate families and bifurcations of periodic orbits, start- 
ing with the study of models based on perturbed oscillators (again, for 
a review, Contopoulos, 2002) and gradually exploring more realistic 
galactic potentials with numerical (Miralda-Escude and Schwarzschild, 
1989; Fridman and Merritt, 1997) and semi-analytical (de Zeeuw and 
Merritt, 1983; Scuflaire, 1995) approaches. 

One of the most powerful analytic tools is the normal form approxi- 
mation of a non integrable system. Although the normal form approach 
is quite widespread in galactic dynamics, its use in studying stability 
of periodic orbits has not been as systematic as the theory could allow 
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(Sanders and Verhulst, 1985). Aim of the present paper is to apply the 
Lie transform normalization method (Dragt and Finn, 1976; Finn, 1984; 
Koseleff, 1994) to approximate the dynamics of a Binney logarithmic 
potential (Binney and Tremaine, 1987). We compare the findings to 
that of Miralda-Escude and Schwarzschild (1989), who employ purely 
numerical techniques to implement the Floquet method and to that of 
Scuflaire (1995), who studies the stability of axial orbits by solving the 
Hill-like perturbation equation with the Lindstedt-Poincare approach. 
Another example that is briefly treated is provided by the galactic 
Schwarzschild (1979) potential with a comparison to the results of de 
Zeeuw and Merritt (1983). These authors based their approach on the 
averaging procedure of normalization: it is therefore interesting a com- 
parison with that method also. We remark that, in careful numerical 
computations, the accuracy of predictions is usually much higher than 
in approximate analytical approaches. However, a reliable analytic tool 
is of invaluable help to gather a global overview of the behavior of the 
system. 

To study the linear orbital stability of the main periodic orbits with 
a truncated normal form one can proceed in essentially two ways: the 
most general and exhaustive is that of determining the explicit form of 
the normal modes and solve the equation of their perturbation. A less 
general but easier approach is that of determining the nature of the 
fixed points on a surface of section. This is constructed with the aid of 
the approximate integral of motion provided by the normalization. The 
first method is in general quite cumbersome and can be applied when 
the procedure of reduction to a single degree of freedom Hamiltonian 
system and the use of action-angle variables lead to a reasonably simple 
system of equations. The second one is clearly less general but relies 
on simple geometric arguments related to the Hessian of a polynomial 
in its critical points and is, at least in principle, quite easy to imple- 
ment. In this work we are going to apply both methods to perform the 
comparison mentioned in the paragraph above. 

In galactic dynamics, the periodic orbits along the axes of symmetry 
(axial orbits) play a particularly important role; moreover they are 
easily identified both as normal modes of the reduced system and as 
"central" fixed points on the surfaces of section. Therefore, we will limit 
the detailed evaluation of the stability characteristics in the parameter 
space to these axial orbits. However, both procedures we have followed 
are quite general and can be directly applied to all periodic orbits of 
sufficiently low commensurability. 

From the results obtained, we can state that the predictive power 
of the normal form ranges well outside the neighborhood in which the 
expansion of the original Hamiltonian is performed. It is rather related 
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to the extent of the asymptotic convergence radius of the approximate 
integrals of motion. However, in concrete applications, the validity of 
the prediction has to be corroborated with an independent evaluation 
of the best suited resonant normal form of the problem at hand. A 
criterion for the choice is illustrated in the last section devoted to the 
applications and is connected to the ratio between the frequency of the 
periodic orbit and that of a normal perturbation to it. 

The plan of the paper is as follows: in Section 2 we recall the proce- 
dure of normalization as applied to reflection symmetry potentials; in 
Sections 3 and 4 we study the 1:1 and 1:2 resonances respectively; in 
Section 5 we reconstruct approximate integrals of motion; in Section 6 
we compare our analytical results with those available in the literature. 



2. The procedure 

2.1. General 

Suppose the original system is given by a Hamiltonian 

H(p,ci) = l(pl+P 2 y ) + V(x 2 ,y 2 ), (1) 

with V a smooth potential with an absolute minimum and reflec- 
tion symmetry with respect to both axes. In our case we will use the 
Schwarzschild and logarithmic potentials described below. We expand 
the potential up to some given degree so that 

oo 

V(x,y;e) = J2^Vn(x,y) (2) 

n=0 

and look for a new Hamiltonian given by 

oo 

K(P, Q;s) = J2 z n K n (P, Q; e) = M^H(p, q; e) , (3) 

n=0 

where P, Q result from the canonical transformation 

(P,Q)=M,(p,q). (4) 
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By expanding (3) in power series of e and equating the coefficients of 
the same order, one has 

Kq= Hq, 



K n — H n + X)m=o M n - m H m — H n + M n H + ^m=i M n ^ m H, 



(5) 

The linear differential operator M ff is defined by 

M g = e- £L ^ e- £2L ^ • • • e~ £ " L ^ • • • , (6) 

where the functions g n are the coefficients in the expansion of the gener- 
ating function of the canonical transformation and the linear differential 
operator Lg is defined through the Poisson bracket 

w={s , „ = _«*£/). 

fr[\dqidpi dpidqj 

The exponentials in the definition of M g are intended as the formal 
sum of a power series so that it gives rise to a near identity coordinate 
transformation known as Lie series. 

The unperturbed part of the Hamiltonian, Hq, determines the form 
of the transformation. In fact, the new Hamiltonian K is said to be in 
normal form if 

{H ,K} = 0. (8) 

This condition is used at each step of the procedure to determine each 
function g n in order to eliminate as much as possible terms in the 
new Hamiltonian. The only terms of which K is made of are those 
staying in the kernel of the operator Lh associated to Hq through the 
definition above. The procedure is stopped at some "optimal" order 
and therefore in all ensuing discussion we refer to a "truncated" normal 
form. Hq must be considered a function of the new coordinates at each 
step in the process: it is therefore an integral of the motion for the new 
Hamiltonian K. The function 

I = K - Hq (9) 

can be therefore used as a second integral of the motion conveying 
approximate informations on the dynamics of the original system. For 
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practical applications (for example to compare results with numerical 
computations) it is useful to express approximating functions in the 
original physical coordinates. Inverting the coordinate transformation, 
the new integral of motion can be expressed in terms of the original 
variables. Denoting it as the power series 

oo 

i = ( 10 ) 

n=0 

its terms can be recovered by means of 

n— 1 

/„ = H n - K n + M n~m[H m - I m ] , n>l. (11) 

m=l 

We remark that in all subsequent applications involving series expan- 
sions, the role of the perturbation parameter can also be played by 
the size of the neighbourhood of the origin where the Hamiltonian 
is considered. Therefore the powers of the parameter e are left in all 
expansion formulas just to indicate their order and are treated as unity 
in the computations. 

2.2. Galactic potentials 

The model potentials we will consider are the Binney logarithmic po- 
tential (Binney and Tremaine, 1987) and the Schwarzschild (1979) 
potential. In both cases we will actually need an expansion of the form 
(2) and we will assume that each term can be written as a homogeneous 
polynomial of degree k + 2 of the form 

V k (x,y) = E^(^)^ HH . (12) 

The logarithmic potential is 

V =Uog{l + x 2 + y 2 /q 2 ) (13) 

and plays a very important role in galactic dynamics because, despite 
its very simple form, it has realistic features like a density distribution 
compatible at large radii with flat rotation curves. The form written 
here is simplified by the choice of fixing the length scale (the "core 
radius" R c ) equal to one, but this is not a limitation due to the invari- 
ance in both the length scale and the energy scale. With these units, 
the energy E may take any non negative value 

< E < oo. (14) 
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The parameter giving the "ellipticity" of the figure ranges in the inter- 
val 

0.6<g<l. (15) 

Lower values of q can in principle be considered but correspond to a non 
physical density distribution. Values greater than unity are included in 
the treatment by reversing the role of the coordinate axes. The series 
expansion of the logarithmic potential is 

v = T, { -^— {* 2 + y 2 iey, (i6) 
3=1 3 

so that the lowest order coefficients are 

a 20 = w\ = !> °02 = u\ = l/q 2 , (17) 
a 40 = -1, a 22 = -l/q 2 , a 04 = -l/<? 4 , (18) 
a 60 = 1, a 42 = 3/g 2 , a 24 = 3/q A , a 6 = l/q 6 - (19) 

The Schwarzschild (1979) potential is to be considered more for its his- 
torical role rather than for its practical usefulness. However, it has been 
deeply investigated and is therefore a good benchmark for comparison. 
It can be written as 

2 2 

V = u{r) + X ~ V w(r) + 1, (20) 



where 



U = ~tog(>/rT^ + r) -ci 2(1+ r c 2 2r2) 3/ 2 > (21) 



„2 



w = - C3 (iw)^' (22) 



are two functions of r = \J x 2 + y 2 such that < u, w < 1 and the c's 
are fixed constants. With the choice of de Zeeuw and Merritt (1983) 1 

ci = 0.064, c 2 = 0.655, c 3 = 0.015, c 4 = 0.481, (23) 

the lowest order coefficients are 

wi = 0.421, w 2 = 0.601, (24) 
a 40 = -0.042, a 22 = -0.174, a 04 = -0.307, (25) 
a 60 = -0.006, a 42 = 0.221, a 24 = 0.460, a m = 0.233. (26) 



1 Note the correction in C4 with respect to the value reported in the appendix of 
de Zeeuw and Merritt (1983), necessary to comply with the other reported constants. 
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The energy range in the Schwarzschild potential is 

< E < 1. (27) 

2.3. Detuning the normal form 

The natural setting in which one can perform a low order normalization 
is therefore that of a perturbed quadratic Hamiltonian with a potential 
starting with a harmonic term. In the general case in which the frequen- 
cies are rationally independent, the kernel of the operator associated 
to 

H = l -{pl+p 2 y +ulx 2 + uly 2 ) (28) 

is trivial, consisting only of functions of the partial energies in the 
harmonic potential: it is customary to refer to the normal form con- 
structed in this "Birkhoff" normal form (Birkhoff, 1927). The 
presence of terms with small denominators in the expansion, forbids 
in general its convergence. It is therefore more effective to work since 
the start with a "resonant" normal form, with a "richer" kernel that 
allows to reconstruct the main natural resonances shaping the phase- 
space of the system. To catch the main features of the orbital structure, 
we therefore approximate the frequencies with an integer ratio plus a 
small "detuning" that we assume 0(e 2 ) 

^ = ™+eH (29) 
uj 2 n 

and we speak of a detuned (m:n) resonance, with m + n the order of 
the resonance. 

We have to put the system in a form suitable to apply the normal- 
ization procedure: we rescale variables in order to put the Hamiltonian 
in the form 

-i oo fc+2 i 

H = l[(m + n6)(plW)+n(plW)]+j:E ^f^V^ (30) 

where we have used the same notation for the rescaled variables and 

, _ "Q(j,fc+2-j) ( , 

°(j,k+2-j) ~ JJ2 (fc+ 4 -j)/2- ^ 

The procedure is now that of an ordinary resonant "Birkhoff-Gustavson" 
normalization (Gustavson, 1966; Moser, 1968) with two variants: the 
coordinate transformations are performed through the Lie series and 
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the detuning quadratic term is treated as a term of higher order and 
put in the perturbation. 

2.4. Choice of the resonance 

Given an arbitrary pair of unperturbed frequencies, it could seem better 
to approximate their ratio as close as possible with a suitable pair of 
integers. However, beside possible computational problems, there are 
arguments on which a more effective choice can be based. Actually, 
the resonance should be of the lowest possible commensurability giving 
rise to a frequency ratio compatible with the dynamics of the actual 
system. The reason for this is that, the lower the order of the resonance, 
the richer the family of terms compatible with it that are available to 
construct the normal form. 

Moreover, another argument in favor of low order resonances comes 
from their role in the stability properties of periodic orbits. A typical 
situation is that in which a family of periodic orbits becomes unstable 
when a low order resonance occurs between its fundamental frequency 
and that of a normal perturbation: the simplest case is given by an 
axial orbit that, depending on the specific form of the potential, can be 
unstable through bifurcation of loop orbits (1:1 resonance), "banana" 
orbits (1:2 resonance), "fish" orbits (2:3 resonance), etcetera. Therefore 
a detuned low-order resonant normal form can be quite accurate in 
describing the corresponding bifurcations. 

Finally, it must be emphasized that the structure of a resonant nor- 
mal form is also affected by the symmetries of the original system. The 
normal form must preserve these symmetries and this in general also 
leads to a criterion for truncation. In the present instance of a double 
reflection symmetry, given a resonance ratio m/n, the normal form 
must contain at least terms of degree 2(m + n) (see, e.g. Tuwankotta 
and Verhulst, 2000). Therefore, the criterion we have adopted in this 
paper has been that of working with the lowest order truncated normal 
form incorporating the symmetries of a typical galactic potential: the 
1:1 symmetric resonance which allows to truncate the normal form to 
degree 4 and the 1:2 symmetric resonance which requires to truncate 
the normal form to degree 6. A systematic investigation of the optimal 
order of truncation has recently been performed by Contopoulos et 
al. (2003) and Efthymiopoulos et al. (2004). Their results confirm the 
rapid decrease of the optimal order with the radius of the phase-space 
domain in which expansions are computed: we may conjecture that if 
we are interested in the global dynamics and accept a moderate level 
of accuracy, with this very conservative approach we can get reliable 
information up to the breakdown of the regular dynamics. 
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3. 1:1 symmetric resonance and first order normalization 

A Lie transform normalization truncated to the second order gives the 
following expression of the first-order normal form (cfr. Belmonte et al. 
2006) 

= ¥(Pi +QD + T2 [MA 2 + Qi 2 ) 2 + MP 2 2 + Q2 2 ) 2 

+ b -§ (P 2 2 (3Pi 2 + Q1 2 ) + Q 2 2 (Pi 2 + 3Qi 2 ) + 4P 1 P 2 Q 1 Q 2 ) , 



K (m:n) = n §{p 2 + Q?) + ^ ^(p^ + Q 2) 2 + ^(p^ + 



2\2 



+fc 2 6 22 (P 1 2 + Q 1 2 )(P 2 2 + Q 2 2 ) 



where fci, k 2 are rational numbers dependent on m and n, Eq. (29) has 
been used, that in the present instance reads 



ui = (l + e 2 8)oj2 



(32) 



and the canonical variables P, Q are as in Eq. (3). We see that in this 
case we have the same situation as in the first order averaging approach 
(Verhulst, 1996): the 1:1 resonance or all other resonances. 

To show more clearly how the symmetries influence the structure of 
the normal form, we have that it can be written as 

K=J 1 + J 2 + £ 2 [5J 1 + - (b m J\ + &04 Jl + |&22 Ji Ji (2 + COS (2#i - 2# 2 ) ) )] 

o 3 

(33) 

where an overall rescaling by a factor u; 2 has been performed and the 
action-angle variables are introduced according to 



Qi = \JlJ\ sin 6*i, 

Pi = V2Jicos6»i, 

Q 2 = V^^sin^, 

P 2 = V2J2COS02. 



(34) 
(35) 
(36) 
(37) 



In fact, inverting these expressions and putting them into (33) we get 

(38) 



lj 2 K = H Q + u 2 e 2 K^ 1) . 



The structure of (33) displays the effect of the symmetries on the res- 
onant part: angles appear only through the combination 29\ — 26> 2 and 
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this shows why the symmetric 1:1 resonance can also be dubbed a "2:2" 
resonance. 

We can use (33) to identify the main periodic orbits. The procedure 
is the following (Sanders and Verhulst, 1985, sect. 7.4). We perform the 
following canonical transformation to "adapted resonance coordinates" 

V> = 2(0! - 2 ), (39) 

X = 2(0 1 + 2 ), (40) 

Ji = (£ + K)/2, (41) 

J 2 = (£-1l)/2. (42) 

In this way we get 

K = -S{£ + K) + A(£ 2 + TZ 2 ) + B£TZ + C{£ 2 - TZ 2 ){2 + cos i>) (43) 

where the new action £ is the additional integral of motion and has 
therefore been subtracted to get the effective Hamiltonian 

K-^. (44) 



The coefficients 



A = J^ 40 + 6 ° 4 )' ( 45 ) 
B = ^(&40-M, (46) 
C = ^b 22 , (47) 

appearing in K, give it the simplest form. 

Considering the dynamics at a fixed value of £ , we have that K 
defines a one-degree of freedom (ip,TZ) system. We get the following 
equations of motion 

ij, = K n = -5 + B£ + 2(A-C(2 + cos4>))1l, (48) 

U = -K f = c(£ 2 - U 2 ^ sin V- (49) 

Let us determine the fixed points of this system: these in turn give the 
periodic orbits of the original system. The right hand of (49) vanishes 
either for 1Z = ±£ or for ip = 0, ±7r. In the first case, the right hand of 
(48) vanishes when 

5 + 2[B ±2(A-C(2 + cosip))}£ = (50) 



celmecIV.tex; 8/02/2008; 20:33; p. 10 



Stability of axial orbits in galactic potentials 



11 



and the two periodic orbits 



■R = £, J 2 = 0, (Type la), (51) 
n = -£, Ji = 0, (Type lb), (52) 



ensue. In the second case, the right hand of (48) vanishes either when 

5 + 2B£ 
4(3C - A) 



Lt W£ A ^ W = 0), (53) 



or when 

K = W^y < 54 > 

The fixed point in (53) determines the "inclined" orbit 

5 + 2(B + 2(3C-A))£ 
Ji = ^c^A) ' (Type n) " (55) 

Note that 

< Ji < £ (56) 

and this range determines the condition for existence of the orbit of 
Type II. The fixed point in (54) determines the elliptic orbit 

The range (56) still determines the condition for existence of the orbit 
of Type III. 

Let us now consider the question of the stability of the periodic 
orbits. In particular, we are interested in what happens in the case 
of axial orbits of Type I: unfortunately, action-angle variables have 
singularities on these orbits and these affect also the adapted reso- 
nance coordinates. However, the remedy is quite straightforward: to 
use a mixed combination of action-angle variables on the normal mode 
and Cartesian variables for the other degree of freedom. The ensuing 
procedure is then first to determine the condition for the normal mode 
to be a critical curve of the Hamiltonian in these coordinates. Second, 
to assess its nature (Kummer, 1977; Contopoulos, 1978; Sanders and 
Verhulst, 1985, sect. 7.4.4): the condition is found by considering the 
function 

if W = K + fiH , (58) 

where \x has to be considered as a Lagrange multiplier to take into 
account that there is the constraint Hq = £ . The Lagrange multiplier 
is found by imposing 

dK^ = (59) 
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on the normal mode. Its nature is assessed by computing the matrix of 
the second derivatives of K^: if the Hessian determinant of the second 
variation is positive definite the mode is elliptic stable; if it is negative 
definite the mode is hyperbolic unstable. 

In the case of the y-axis orbit of Eq. (52), good coordinates are 

Qi = X, (60) 

Pi = U, (61) 

Q 2 = V2Jsm9, (62) 

P 2 = v^cosfl, (63) 

so that the periodic orbit is given by 

X = U = 0, J = 8. (64) 
The terms in the normal form are then 

H = ^(X 2 + U 2 ) + J (65) 

and 

K 2 = \5{X 2 + U 2 ) + l 2 [b m (X 2 + U 2 ) 2 + 4&04 J 2 ] + 

±b 22 J [2(X 2 + U 2 ) + (X 2 - U 2 ) cos 29 + 2XU sin 26} . [ ' 

It is straightforward to check that, in this case, imposing Eq. (59) on 
the periodic orbit defined by Eq. (64), we get the equation 

^+1 + ^04^ = 0, (67) 

which allows to find the required value of the Lagrange multiplier. With 
this result, the matrix of the second derivatives of on the normal 
bundle to the periodic orbit is 

If 88 + 8 [b 22 (2 + cos 20) - Qb m ] 8b 22 sin 29 

8\ 8b 22 sm26 85 + 8 \b 22 (2 - cos 26) - 66 04 ] 

(68) 

The equation detif w (8) = gives 

(36^ 4 - 24604^22 + 3b 2 2 )8 2 - 32(36 04 - b 22 )58 + 645 2 = (69) 
with roots 

85 ^ 85 
66 04 - »22 3(26 04 - 022) 

Recalling the rescaling in (33), the physical energy is given by 

E = uj 2 8. (71) 
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If, as in the application cases that will be examined later, the first 
coefficient in Eq. (69) is positive, the range of instability of the y-axis 
orbit is 

8uj 2 5 „ 8tJ2<5 , , 

<E< — 2 , (72) 



66 04 - b 2 2 " " 3(26 4 - b 2 2) ' 

Proceeding in the same way in the case of the x-oxis orbit of Eq. (51), 
the analogous expressions 

g ^ <E< (73) 
3 b 2 2 - 26 40 622 - 6640 

can be written for the other conditions. 



4. 1:2 symmetric resonance and second order normalization 

In the case of the 1:2 resonance in presence of reflection symmetries 
about both axes, we know that the normal form must be pushed at 
least to degree 2 x (1 + 2) = 6. We therefore have to perform a further 
step of normalization and include K4 in the normal form. Its expression 
is quite involved (cfr. Belmonte et al. 2006), but we can exploit the 
change of variables to action-angle coordinates to see that the normal 
form has the structure 

K = J 1 +2J 2 +e 2 (25J 1 +P 2 (Ji, J 2 ))+e 4 (P 3 (Ji, J 2 )+kJlJ 2 cos(40i-20 2 )), 

(74) 

where Eq. (29) has been used, that in the present instance reads 

^i=Q + e 2 ^w 2 , (75) 

the polynomials P 2 and P3 are homogeneous of degree 2 and 3 respec- 
tively 

3 1 3 

Pi = g&to^i + -gb 22 JiJ 2 + -6 4«/| 

P3 = -3^ ((1026l - 1606 60 ) Jf + (106 2 2 + 72&22&40 - 96b 42 )J 2 J 2 
+(36604^22 + 106 2 2 - 966 24 ) JiJl + (516g 4 - 1606 06 ) J 2 3 
and 

k = -^92 (3^2-3622640-8642). (76) 
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The study of existence and stability of periodic orbits proceeds in the 
same way as in the previous section: we limit the presentation to the 
axial orbits. In the present case they are given by 

n = 28, J 2 = 0, (Type la), (77) 
K = -8/2, Ji = 0, (Type lb). (78) 

The procedure to determine the condition of stability of axial orbits 
lead to analyze critical curves of the function (58) where K is given by 
(74). Considering the x-axis (type la, Eq. (77)) orbit, good coordinates 
are 

Qi = v^sinfl, (79) 
Pi = VzJcosO, (80) 
Q2 = Y, (81) 
P2 = V, (82) 

so that the periodic orbit is given by 

Y = V = 0, J = 8, (83) 

and 

H = J + Y 2 + V 2 . (84) 
The condition Eq. (59) allows to find the Lagrange multiplier 



H = -1-2S-- 
4 



(85) 



The equation detK^ (8) = 0, obtained by computing the matrix of 
the second derivatives of on the normal bundle to the periodic 
orbit (83), is 



{A 2 8 2 + A X S + 7685)(A' 2 8 2 + A X S + 7685) = 0, (86) 



where 



A\ = 48(66 40 - 622), (87) 
A 2 = 2b 2 22 + 33& 2 2&40 - 306&1 - 566 42 + 4806 60 , (88) 
A' 2 = 2b 2 2 + 396 22 6 4 o- 3066l - 406 42 + 4806 60 . (89) 

We see that the inequalities to be satisfied in this case are quadratic in 
the energy, contrary to the linear case of (72) and (73). Below we discuss 
the specific example of the x-axis orbit in the logarithmic potential. 
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5. Approximate integral in the original variables 

The general procedure illustrated so far can be applied to study the 
dynamics associated to every reduced normal form. Specifically, it can 
be extended to forms truncated to an arbitrary high order and to all 
normal modes and periodic orbits associated to the given resonance. 
However, if the 1:1 and 1:2 cases treated here are rather simple, com- 
putations become quite cumbersome for higher orders. It can be useful 
to exploit an alternative, less general, but easier procedure. Recalling 
the generic expression of the second invariant (11), if we truncate at 
first order, we have 

h = H + e 2 (V 2 - K 2 ). (90) 

This is the best approximate integral of motion of a symmetric per- 
turbed 1:1 oscillator to order e 4 in the perturbation parameter: in fact, 
due to the symmetries of the problem, odd-degree terms are absent 
both in the normal form and in the approximate integral. To assess the 
stability of, e.g., the y-axis orbit, we may proceed in this way: to use 
I 2 and the conserved energy 

E = \{pl + p 2 y ) + V (x 2 ,y 2 ) + e 2 V 2 (x 2 ,y 2 ) (91) 

to construct an x — p x Poincare section by means of the intersection 
of the function l2(x,y,p x ;E) with the y = hyperplane and the level 
curves of the function 

F = I 2 (x,0,p x ;E); (92) 

to study the nature of the origin as a critical point determining if it is 
either an extremum (elliptic fixed point = stable periodic orbit) or a 
saddle (hyperbolic fixed point = unstable periodic orbit) by using the 
Hessian determinant 

Ep x p x F X x ~ E xpx - (93) 

Clearly, for resonant periodic orbits even the location of the critical 
points can be already quite difficult and this limits the generality of 
the approach. However, in the case of axial orbits, the approach is 
straightforward: in the case of the y-axis orbit, we have that the second 
derivatives in the origin are 

E PxPx = 4u, 2 <5 - 3E ((1 + 5)b 04 - h 22 ^j , (94) 

F xx = 4lo 2 5 - E f 3(1 + 5)5o4 - ^22) • (95) 
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The mixed second derivative F xpx is identically zero in the origin. The 
range of instability (namely, the range where the two second derivatives 
have different sign) is 

^ < /••• - ~, l 06) 



6604(1 + 5) - b 22 ' 3(26 04 (1 + 5) - b 22 ) 

Comparing with (72) we have an order 5 2 discrepancy whose origin is 
connected with the construction of the function (92) where the term 
with the detuning is not considered as a perturbation. Analogous in- 
equalities can be written for the x-axis orbit by constructing a y — p y 
Poincare section and studying the level curves of the function 

F = I 2 (0,y,p y ;E) (97) 

obtained by means of the intersection of the function I 2 (x,y,p y ; E) 
with the x = hyperplane. 

On the other side, we can develop this parallel approach of deter- 
mining the nature of fixed points on the surface of section for higher 
order approximate integrals too. Using (11), if we truncate at second 
order, we have 

h = H + e 2 {V 2 - K 2 ) + e\V 4 - {g 2 , K 2 } - K A ) (98) 

where g 2 is the first order generating function determined in the first 
step of the normalization. This is the best approximate integral of 
motion of a symmetric perturbed 1:2 oscillator to order e 6 in the per- 
turbation parameter. Let us consider the x-axis orbit: the construction 
of the y — p y Poincare section and the study of the critical points of 
the function 

F = h(0,y, Py ;E), (99) 

obtained by means of the intersection of the function I^x, y,p y ; E) with 
the x = hyperplane, proceeds as above. Concerning the nature of the 
fixed point in the origin (linked to the stability of the x-axis orbit), we 
get quadratic inequalities in the energy analogous to those arising from 
(86): the discrepancy between the two methods is still of order 5 2 and 
is discussed in the subsequent section about a specific example. 



6. Comparison of analytical and numerical results 

We have applied the theory discussed in the previous sections to the 
logarithmic and the Schwarzschild potentials. We can compare these 
analytical predictions concerning the stability energy ranges of axial 
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orbits with other results, mainly numerical, from the literature. In par- 
ticular, we have chosen the work by Miralda-Escude and Schwarzschild 
(1989, MES in what follows) who made an accurate numerical explo- 
ration of the phase-space structure of the logarithmic potential: they 
give the existence parameter ranges of the main families of periodic 
orbits and the bifurcation ensuing from instability thresholds in two 
models (corresponding to ellipticity values q = 0.7 and q = 0.9) deter- 
mined by solving the perturbation equations with the Floquet method. 
The authors use the core radius R c to parametrize the sequence of 
periodic orbits because they were interested in comparing the results 
with the case of the singular scale-free case R c = 0: in fact they fix 
the energy E = in all their computations and vary R c and q. For 
our purposes, it is more natural to use the energy as a parameter and 
therefore we have made the conversion 

E = -\ogR c (100) 

to compare the results. Another approach to the study of the stability of 
axial orbits in the logarithmic potential has been followed by Scuflaire 
(1995, SC in what follows) who solved the Hill-like equation of the 
perturbation by means of a Poincare-Lindstedt series expansion up to 
order 20. SC uses a, the amplitude of the axial orbit, as a parameter: 
the conversion to the energy is given by 

E = Uog{l + a 2 ). (101) 

In tables I. and II. we display the results in the case of the logarithmic 
potential: in the first column are the two values of the parameter q for 
which reliable numerical estimates are available; in the second are the 
analytical predictions according to the procedures of Section 3 for table 
I and 4 for table II; in the third are the analytical predictions according 
to the procedures of Section 5; in the fourth the numerical prediction 
of MES; in the last the semi-analytical prediction of SC. 

The values shown in table I. represent the energy thresholds to 
instability of the short y-axis orbit: this orbit becomes unstable when 
its frequency falls to that of a normal perturbation and as soon as this 
happens a loop orbit bifurcates. The values of the second and third 
columns are therefore computed using the results of the 1:1 resonance 
treatment of Sections 3 and 5, namely the first inequality in (72) and the 
first in (96). The computations are performed by using the coefficients 
given in (18) recalling definition (31). Clearly the accuracy of both 
predictions is better for q = 0.9 corresponding to a smaller detuning of 
the 1:1 resonance. Overall, we see that the prediction using the invariant 
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Table I. Energy threshold of instability for the short 
y-axis orbit in the logarithmic potential 



9 


Normal — Form 


Invariant 


MES 


SC 


0.7 


0.52 


0.86 


0.72 


0.80 


0.9 


0.19 


0.22 


0.21 


0.25 



Table II. Energy intervals of instability for the long x-axis orbit in 
the logarithmic potential 



g 


Normal — Form 


Invariant 


MES 


SC 


0.7 
0.9 


1.32 ^4.82 
0.85 


1.71 ^4.86 
0.80 


1.52-^4.29 
3.64 


1.42 ^ - 



curves performs significantly better than that obtained directly working 
with the normal form. 

The long x-axis orbit can suffer instability, but not through a 1:1 
resonance. This orbit becomes unstable when its frequency falls to 1/2 
of that of a normal perturbation and, as soon as this happens, a banana 
orbit bifurcates. Therefore, we turn the attention to table II. where the 
values of the second and third columns are computed using the results 
of the 1:2 resonance treatment of Sections 4 and 5: namely, the con- 
dition in order to the quartic in (86) is negative, using the coefficients 
computed from (18) and (19) and the nature of the fixed point in the 
origin as determined through function (99). An important prediction 
of the numerical approaches is that, if the ellipticity is less than a given 
value, there is an instability interval in energy: actually in the figure of 
SC it is not possible to determine the upper extremum of it, but this 
interval is very well determined in the computations of MES. We see 
that the analytical estimates of the instability interval are not very far 
from them, with an acceptable prediction of the upper boundary of the 
interval too. Our treatment also provides a critical value of q c = 0.78 



Table III. Energy threshold of instability for the short y-axis orbit in 
the Schwarzschild potential 



Normal — Form 


Invariant 


DZM averaging 


DZM numerical 


0.19 


0.28 


0.21 


0.27 
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of the ellipticity, beyond which the interval disappears. It would be 
interesting to numerically test this prediction. For higher values of q the 
analytical estimates of the threshold of instability are not so good and 
the reason is that the detuning of the 1:2 resonance becomes excessively 
large. 

Finally, in table III. we report the comparisons concerning the short 
y-axis orbit in the Schwarzschild potential, with the coefficients listed in 
(24) and (25): in the third column we report the prediction given by de 
Zeeuw and Merritt (1983, DZM) using the first order averaging of the 
equations of motion and in the fourth the corresponding prediction of 
the same authors using a numerical integration. We see that, as usual, 
the analytical prediction based on the invariant is quite accurate. 

We conclude observing that applications like those presented here 
confirm that this very conservative strategy provides sufficient quali- 
tative and quantitative agreement with other more accurate but less 
general approaches. An important issue to clarify is if, in order to im- 
prove the degree of accuracy, it is enough to simply truncate the normal 
form to higher orders or if it is rather necessary to incorporate in a more 
careful way the detuning in higher order terms of the perturbation. 

Acknowledgements 

Our thanks to Giuseppe Gaeta and Ferdinand Verhulst for very fruitful 
discussions. 

References 

Belmonte C. Boccaletti D. and Pucacco G. (2006) "Approximate First Integrals for 
a Model of Galactic Potential with the Method of Lie Transform Normalization" , 
submitted to the Proceedings of the Saarifest 2005, E. Perez-Chavela and J. Xia 
editors. 

Birkhoff G. D. (1927) Dynamical Systems, Amer. Math. Soc. Coll. Publ., 9, New 
York, USA 

Binney J. and Tremaine S. (1987) Galactic Dynamics, Princeton University Press. 
Boccaletti D. and Pucacco G. (1999) Theory of Orbits: Vol. 2, Springer- Verlag, 
Berlin. 

Contopoulos G. (1978) "Higher order resonances in dynamical systems" . Cel. Mech., 
18, 195-204. 

Contopoulos G. (2002) Order and Chaos in Dynamical Astronomy, Springer- Verlag, 
Berlin. 

Contopoulos G. Efthymiopoulos C. and Giorgilli A. (2003) "Nonconvergence of 
formal integrals of motion". J. Phys. A: Math. Gen., 36, 8639-8660. 

de Zeeuw T. and Merritt D. (1983) "Stellar orbits in a triaxial galaxy. I. Orbits in 
the plane of rotation". Astrophys. J., 267, 571-595. 



celmecIV.tex; 8/02/2008; 20:33; p. 19 



20 



C. Belmonte, D. Boccalctti and G. Pucacco 



Dragt A. and Finn J. M. (1976) "Lie series and invariant functions for analytic 

symplectic maps". J. Mat. Phys., 17, 2215-2227. 
Efthymiopoulos C. Giorgilli A. and Contopoulos G. (2004) "Nonconvergence of 

formal integrals: II. Improved estimates for the optimal order of truncation". 

J. Phys. A: Math. Gen., 37, 10831-10858. 
Finn J. M. (1984) "Lie series: a perspective. Local and Global Methods of Nonlinear 

Dynamics". Lecture Notes in Physics, 252, 63-86. 
Fridman T. and Merritt D. (1997) "Periodic orbits in triaxial galaxies with weak 

cusps". Astron. J., 114, 1479-1487. 
Gustavson F. (1966) "On constructing formal integrals of a Hamiltonian system 

near an equilibrium point". Astron. J., 71, 670-686. 
Koseleff P. V. (1994) "Comparison between Deprit and Dragt-Finn perturbation 

methods". Cel. Mech. & Dynam. Astron., 58, 17-36. 
Kummer M. (1977) "On resonant Hamiltonians with two degrees of freedom near 

an equilibrium point". Lecture Notes in Physics, 93, 57-75. 
Miralda-Escude J. and Schwarzschild M. (1989) "On the orbit structure of the 

logarithmic potential". Astrophys. J., 339, 752-762. 
Moser J. (1968) "Lectures on Hamiltonian systems". Mem. Am. Math. Soc, 81, 

1-60. 

Sanders J. A. and Verhulst F. (1985) Averaging Methods in Nonlinear Dynamical 

Systems, Springer- Verlag, New York. 
Schwarzschild M. (1979) "A numerical model for a triaxial stellar system in 

dynamical equilibrium" . Astrophys. J., 232, 236-247. 
Scuflairc R. (1995) "Stability of axial orbits in analytic galactic potentials". Cel. 

Mech. & Dynam. Astron., 61, 261-285. 
Tuwankotta J. M. and Verhulst F. (2000) "Symmetry and Resonance in Hamiltonian 

Systems". SI AM J. Appl. Math., 61, 1369-1385. 
Verhulst F. (1996) Nonlinear Differential Equations and Dynamical Systems, 

Springer- Verlag, Berlin. 



celmecIV.tex; 8/02/2008; 20:33; p. 20 



